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We apply numerical simulations to study the criticality of the 3D Ising model 
with random site quenched dilution. The emphasis is given to the issues not 
being discussed in detail before. In particular, we attempt a comparison of 
different Monte Carlo techniques, discussing regions of their applicability and 
advantages/disadvantages depending on the aim of a particular simulation set. 
Moreover, besides evaluation of the critical indices we estimate the universal 
ratio r+/r^ for the magnetic susceptibility critical amplitudes. Our estimate 
r^/F" = 1.67 ±0.15 is in a good agreement with the recent MC analysis of 
the random-bond Ising model giving further support that both random-site and 
random-bond dilutions lead to the same universality class. 

PACS: 61.43.Bn, 64.60.Fr, 75.10.Hk 

1. Introduction 

Since its introduction in the early 1920-ies the Ising model serves as a paradigm 
to study criticality of interacting many-particle systems where the single particle 
can be in two possible states. Such a general problem statement supported in the 
1960-ies by the universality and scaling hypothesis led to the present situation, 
when the Ising model is used to explain criticality and scaling in such different fields 
as ferromagnetism and binary mixtures from the one side or networks and text 
series from the other side. Subsequently, the random-site Ising model is of primary 
importance to understand the influence of the structural disorder on criticality. And 
this explains a lot of theoretical, experimental and computational efforts invested so 
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far to shed light on the phenomena happening there. 

In particular, the three-dimensional (3D) Ising model subject to a weak dilution 
by non- magnetic impurities changes its universality class [ . Already this phe- 
nomena called for its explanation and verification. We address the reader to recent 
reviews [ |2l E] where theoretical, experimental and computational data is collected 
and compared. In our study, we use the Monte Carlo (MC) approach and address the 
comparison of three different simulational techniques analysing critical behaviour of 
the 3D quenched diluted Ising model. Some of the outcome is obvious, but we think 
it is instructive to perform such simulations now, when certain MC tools are some- 
times discarded for not obvious reasons. Moreover, a part of our simulations concerns 
the critical amplitude ratios of the 3D random-site Ising model - a question which 
was not intensively analysed so far, and where a disagreement between the MC and 
theoretical predictions exists. 

It is our big pleasure and honour to contribute this paper to the Festschrift 
devoted to Reinhard Folk's 60th birthday: his contribution to the theory of phase 
transitions in general and to the problem we consider here is hard to be overesti- 
mated. 

2. Simulation details and results 

The Hamiltonian of the random-site Ising model has the following form 

H = - JY^ CiCjSiSj, (1) 

where (ij) denotes the summation over the nearest neighbour sites of 3D simple 
cubic lattice, J is the interaction constant, q = 1 if the i-th site is occupied by spin 
and Cj = otherwise, the Ising spins Si take two values ±1. Occupied sites (c, = 1) 
are considered to be randomly distributed and quenched in a fixed configuration. 
For every observable, discussed below, first the Boltzmann average with respect to 
the spin subsystem is performed for the fixed site configuration, subsequently the 
averaging over different disorder realisations (configurational average) is performed. 
We will use the following notations, the number of all sites is = L'^ and the 
number of sites carrying a spin is Np. The concentration of spins is defined therefore 
as p = Np/N. 

Several separate sets of simulations are performed in this study with the aim to 
investigate different characteristics of the model. In all the simulations the concentra- 
tion of spins was taken fixed and equal to p = 0.85. For the 3D random site Ising 
model, at spin concentration p ~ 0.8 the correction-to-scaling terms appear to be 
particularly small [|3]. Therefore, for the concentration chosen we do not take these 
terms into account in our further analysis. The simulations were performed on a 
set of following lattice sizes L = 10, 12, 16, 24, 32, 48, 64, 96 with periodic boundary 
conditions. The other details, particularly the number of disorder realisations, the 
simulation lengths and the algorithms employed have been chosen depending on the 
aim of particular simulation set. 
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The aim of the first set of simulations is to estimate the critical temperature of 
the model at different L. Due to the finite-size scaling theory the finite system 
of linear size L will demonstrate an evidence of a critical behaviour at a certain 
temperature Tc{L) which differs from the critical temperature of the infinite system 

T,(L) = T,(oo) + aL-i/-^ + . . . , (2) 

where the correction-to-scaling terms have been omitted. To find Tc{L) the following 
procedure was followed. The initial configuration was prepared by scattering diluted 
sites randomly on the lattice and the original orientations of spins were chosen 
randomly. First, the runs on the smallest lattice size L = 10 were performed on 
10^ disorder realisations during 5-10^ MC sweeps using Metropolis algorithm [|7]. 
The preliminary estimate for the critical temperature was taken from the mean-field 
approximation : 

T^'F = p . T^""^, (3) 

where T^"*"^ is the critical temperature of the pure 3D Ising model [ 6j. During the 
simulations the histograms for the potential energy (0) and for the absolute value 
of the order parameter 

1 ^ 

M = ^\T.c^S,\ (4) 
^^p i=i 

were built and then, using the histogram reweighting technique the temperature 
region around T^^^ was explored. For a given disorder realization, the susceptibility 
was evaluated according to the fluctuational relation 

X = KN,{{M')-{M)'), (5) 

where (. . .) denotes Boltzmann averaging and K = (3 J = J/ksT is the dimensionless 
coupling. The temperatures where x has a maximal value were then averaged over 
all disorder realisations (hereafter, such configurational averaging will be denoted 
by . . .) and used as the working estimate for the critical temperature Tsim{L) for 
L = 10. This temperature was also used as the preliminary estimate for the critical 
temperature for the next lattice size L = 12. Again, for L = 12 we performed 
short runs on 10^ disorder realisations during 5 ■ 10'^ MC sweeps. This provides 
us with the working estimate for the critical temperature for L = 12, Tsjrrt(12). 
The process is continued until all the lattice sizes are processed. The temperatures 
Tsim{L) obtained in this way have been used for the final simulations of 10'^ disorder 
realisations from which all the main results are driven. It is to note here that one 
has different possibilities to define the critical temperature of a finite-size system. 
Along with the definition used in this study, one can define Tc{L) as the temperature 
where the configurational average of x has its maximum. However the same value 
of Tc is expected in the thermodynamic limit and the same scaling behaviour holds 
when approaching this limit. 

The large clusters of correlated spins exist in the vicinity of the critical point, 
that in turn leads to the effect of critical slowing down [ElCSl- As the result, the 
relaxation time of the system increases dramatically. One can estimate the typical 
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relaxation time during the simulation by observing the (configurationally dependent) 
autocorrelation function [ for the potential energy (or, alternatively, for some 
other characteristic of the system) 



where E{t) is the instant value for the energy at some effective time t (in MC 
simulations t is measured in MC sweeps) and to is some time origin. The C(t) is 
averaged over different time origins in a course of simulations. At times large enough 
C(t) decays exponentially according to Debye law 



more details can be found elsewhere [E]. In ((Tj) r has a meaning of the characteristic 
autocorrelation time for a given disorder realisation. Besides the influence on critical 
dynamics (which is beyond the scope of this study), the critical slowing down has 
also some practical implications. Firstly, the estimate of r is vital for determination 
of the required length of the simulation. Only the configurations separated by a 
number of MC sweeps of order r can be considered as statistically independent, and 
the required length of the simulation run should be measured in r scale but not in 
MC sweeps [E]- Secondly, the critical slowing down leads in Metropolis (or other 
local dynamics) algorithms to a dramatic increase of r. This is, of course, due to 
a high energy penalty required to flip a single spin (or pair of spins) in a cluster 
of uniformly oriented neighbours. To overcome this problem, a number of cluster 
algorithms with non-local dynamics have been suggested, with the Swendsen-Wang 
[ IT8] and Wolff [ 14J ones being most commonly used. 

The second set of our simulations is targeting upon the detailed comparison of 
autocorrelation times and of the efficiency of different MC algorithms at a number 
of lattice sizes L for the model under consideration. To this end we performed the 
simulations of 10-15 disorder realizations for each L at the temperatures Tsim{L)- 
We used three different algorithms, the Metropolis, Swendsen-Wang and Wolff ones. 
The simulation length was fixed to 10^ for all cases. For each given disorder realisa- 
tion the value of r has been evaluated according to ((Tj). These values were averaged 
then over all disorder realisations. The results for the average autocorrelation time 
te are presented in tab. 1. The relaxation times are given in MC sweeps, where one 
sweep corresponds to refreshing of each spin. This means, that for the Wolff algo- 
rithm we followed the standard procedure of the autocorrelation time normalisation 
taking into account the average cluster size [El- One can see the different rate of 
pseudo dynamics for all three algorithms used and the autocorrelation time te can 
be seen as a measure of the minimal number of MC sweeps separating two adjacent 
uncorrelated configurations. However, all three algorithms differ significantly on a 
time spent for one MC sweep and we will be more interested in the time required to 
generate the next uncorrelated configuration. This will reflect the true speed of each 
algorithm. The latter can be estimated from an inverse time ttriai spent on Ntriai 



C{t) 



{Eito)E{to + t)) - {E{to)){E{to)) 
{Eito)Eito))-{Eito)){Eito)) 



(6) 



C(t) = exp(-t/r). 



(7) 
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MC sweeps and taking into account that only Ntriai/TE of these are uncorrelated 

^ trial 



L 


10 


12 


16 


24 


32 


48 


64 


96 


te, Metropohs 


5.79 


7.87 


12.9 


26.5 


44.9 


93.1 


160 


336 


te, Swendsen-Wang 


3.55 


3.98 


4.69 


5.54 


6.10 


7.50 


8.55 


10.7 


te, Wolff 


1.18 


1.24 


1.31 


1.46 


1.54 


1.71 


1.81 


1.98 



Table 1. The energy autocorrelation time te for different lattice sizes L measured 
in MC sweeps (see the text for details). 

To this end we performed short trial simulations for 10 disorder realisations each 
of which was well equilibrated, with Ntriai = 10^. The simulations were performed 
for all lattice sizes L = 10 — 96 and using all three MC algorithms. The results are 
presented in fig. One can conclude the following. Both cluster algorithms have 
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Figure 1. Log-log plot for relative speed measured as on inverse time spent on 
generation of one uncorrelated configuration for lattice sizes L = 10 — 96 and 
three different MC algorithms 

significantly higher real speed than the Metropolis one for all lattice sizes. Concern- 
ing the larger lattices examined both cluster algorithms demonstrate comparable 
speed with the Wolff one being slightly faster for L < 80. At L ~ 80 the speed of 
both cluster algorithms equals, and for the larger lattice sizes the Swendsen-Wang 
algorithm looks more preferable (however, we have not done any simulations be- 
yond L = 96 to check this). According to these findings we can conclude that for 
the model and lattice sizes examined in this study the Wolff algorithm is the most 
preferable one and therefore it was used to obtain all the main results in this study. 
We also should note that our realisation of the Metropolis algorithm does not use 
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any special tricks to speed-up the simulation of the Ising model (energy maps, multi- 
spin coding, etc.). Therefore, our estimate of a true speed of Metropolis algorithm 
is subject to the particular code being used. 

The third set of our simulations aimed on the evaluation of the critical indices 
and the extrapolation of the critical temperature Tc(oo) for the infinite system. 
Along with the existing MC estimates of the critical indices of diluted 3D Ising 
model (see e.g. Refs. [013 for reviews) for the sake of consistency we performed our 
independent evaluation of these quantities. To this end the MC simulations of 10^ 
disorder realisations were performed for lattice sizes L = 10, 12, 16, 24, 32, 48, 64, 96 
using Wolff algorithm. The simulation length was chosen equal to 250te for the 
equlibration and IO^te for the production runs. These numbers were chosen to bal- 
ance the errors originated from the limited number of disorder realisations Ndis and 
the finite length Nmc of the simulations, both terms appearing in the total error 
estimate for some quantity A [ IT^ 

AA=(^ + (l + 2r,)^^i^)"\ (9) 



where 



1/2 



<yiA) = {{A^)-{A) (10) 



is the standard deviation of the set of N^is averages {A) calculated within each 
disorder realisations, and 

aA=(m^Wf" (11) 

is the standard deviation of the whole set of NiHsNmc data for the quantity A. The 
ta in Eq.Q is the integrated aucorrelation time r^""* [^] which, in general differs 
from the exponentional one r^^^ introduced in ((7j) in the case of the energy. However, 
for most of realistic models r^"* < r^^^ [ ^] so that the use of r^^^ only which we 
employ in this study is reasonable. 

In evaluation of the critical indices we followed the standard finite-size scaling 
scheme [E]. According to it, the Binder's cumulant 

C/ = l-^ (12) 

can be evaluated as a function of a coupling K for each disorder configuration and 
the maximum value for the slope of this function varies with the system size as L^^'^ 
anywhere in the critical region. In the terms of temperature derivative one has 



dU 

Ik 



dT 



aL'/", (13) 



where the histogram reweighting technique was used to evaluate the cumulant U in 
the neighbourhood of Tc{L) and the numerical derivation was employed. Here and 
thereafter the observables (e. g. ()13p - ()15|) ) are averaged over disorder realisation. 
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The inverse values for v can be found by the hnear interpolation of the log-log plot 
for the values 4^ vs L (see, fig. EI)- The same finite size scahng behaviour is valid 

""^ max 

for the number of logarithmic derivatives for the powers of the magnetisation [[Hj. 
These can be used as additional estimates for the critical index z/, and we evaluated 
two of them, namely 



Ain(M) = -T2^1n(M), (14) 



Ain(M2) = -T2Ain(M2). (15) 



The log-log plots for these derivatives are shown in fig. El too. 




Figure 2. Log-log plots for the maximum values of the configurationally averaged 
derivatives of Binder cumulant (|13|) (discs), ln(M) ()14() (triangles) and In(M^) 
((T3)) (diamonds) 

First of all we would like to stress that all the data points presented in this figure 
interpolate extremely well by the appropriate linear dependencies. This can be seen 
as an other proof that the number of disorder realisations and the simulation length 
that have been used are reasonable. The linear interpolation for the curves presented 
in fig. El leads to the values presented in tab. 2. 

The average of all the three results for the index v shown in tab. 2 yields 

1/ = 0.662 ±0.002 (16) 

Obtained value for v distingtly differs from the theoretical estimate of the corre- 
sponding exponent of the pure 3D Ising model: v = 0.6304(13) [Cni- However it is 
less than the asymptotic value of u for the random 3D Ising model as obtained by 
the field-theoretical RG approach u = 0.678(10) [EHj and by the MC simulations 
ly = 0.6837(53) [m, 1/ = 0.683(3) [HE]. This is an evidence of the fact that for the 
spin concentrations and lattice sizes chosen here the system still crosses over to the 
asymptotic regime. Note, that close estimates for u were found at similar system 
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interpolation of 




V 


dU 

max 


max 
m,ax 


1.507 ±0.008 
1.512 ±0.001 
1.514 ±0.001 


0.664 ±0.004 
0.661 ±0.001 
0.660 ±0.001 



Table 2. The values for the critical index v obtained via linear interpolation of 
the log-log data for the slope of Binder's cumulant and logarithmic derivatives of 
(M) and (M^) 

parameters (and neglecting correction-to-scaling terms) in the MC simulation of the 
random site 3D Ising model: p = 0.9, L = 64 128, v = 0.6644(15) [|3]; as well as of 
the random bond 3D Ising model: ptonds = 0.7, u = 0.660(10) [ITBj. 

The standard finite-size scaling concepts have been also employed for the evalu- 
ation of the other critical indices, particularly these for the susceptibility 7 and for 
the magnetisation f3. Due to this scheme [IH] for finite size system of dimension L in 
the critical region one has 

where F and B are the critical amplitudes and the correction-to-scaling terms have 
been omitted. In (|T7jl Xmax is the averaged maximum value for the susceptibility that 
is achieved at some temperature and (M) is the averaged value of the magnetisation 
at the same temperature. We employed the following scheme. For each disorder 
realisation the temperature T* was found where the x achieves its maximal value 
Xmax- Then the magnetisation (M) is evaluated at the same T*. Afterwards, the 
values of Xmax and (M) have been averaged over all disorder realisations to be used 
for fitting the relations (fTTj) . The results of fittings are presented in figs. I3I4I in a 
form of the log-log plots for the Xmax and for (M) vs L, respectively. 

One should note near perfect linear fits achieved in both cases, leading to the 
following results for 7 and P 



7/1/ = 1.986 ±0.001, (18) 

p/u = 0.509 ±0.001, (19) 

7 = 1.314 ±0.004, (20) 

P = 0.337 ±0.001. (21) 



For the reasons explained above, taking exponent u as an example, the values ()18|)- 
(PTj) slightly differ from their asymptotic counterparts P = 0.3546(28), 7 = 1.342(10) 
[ 0] . They are in a good agreement with the exponents of the 3D Ising model with 
random bond dilution: P/u = 0.515(5), 7/z/ = 1.97(2) [E] giving one more proof 
that both models belong to the same universality class. 

With the value for the critical index u being evaluated in p6j) one can obtain 
the estimate for the critical temperature for the infinite system Tc(oo). It is related 
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Figure 3. Log- log plot for the maximum value of the susceptibility Xmax vs the 
linear size of the system L in the critical region 
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Figure 4. Log-log plot for the magnetization (M) vs the linear size of the system 
L in the critical region 

to the critical temperatures of the finite system of hnear dimension L, Tc{L), via 
the relation Q. The values for Tc{L) were estimated during the same set of simu- 
lations that have been used for the evaluation of critical indices. For each disorder 
realisation we found the sets of temperatures T^{L) where the maximum values for 
the susceptibility (for the i = 1 case) and for the derivatives (fT^ - (fT3j) (the cases 
i = 2,3,4, respectively) are achieved. The data sets T^{L) should be plotted then 
vs the scaled system size L"^/'^ and, ideally, for all i should extrapolate at L — > cxo 
to the same value of Tc(oo). The results are shown in fig. El where we used the data 
for all the lattice sizes from L = 10 to L = 96. As the result of the extrapolation 
procedure we obtain 

r 3.756931 ±0.000058, 
i, ,1 3.756783 ± 0.000243, 
^A^) - S 3.756481 ± 0.000086, 

[ 3.756543 ± 0.000083, 





XI max 

dU / dK\,^nx {oo\ 
d\n{M)/dKU. ^ ' 

d\n(M^)/dK\ 

max 
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Tc{L) 
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3.81 


* ln(M) 
T ln(i\/2) 

A 




3.79 






3.77 


■ 
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lit \ ^ 




3.73 




1 1 



0.00 0.01 0.02 0.03 i\lv 

Figure 5. The linear fits for the critical temperatures TJ:{L) vs the scaled system 
size L~^/'^ for the lattice sizes from L = 10 to L = 96, the index i corresponds to 
the maximum values for different quantities, see notations in the figure 

One should note very good convergence of the extrapolated values for T^*(oo) 
found by examining of different quantities and given in Eg. ()22|1 . The final estimate 
for the critical temperature of the infinite system Tc{oo) can be found as an average 
on all four numbers in Eq. ()22p 

Tc(oo) = 3.7566845 ±0.0001175, /^^ = T-^(oo) = 0.2661922 ± 0.0000083. (23) 

As it is clear already from simple mean- field arguments (see Eq. (3)) our estimate 
for Tc made for the spin concentration p = 0.85 should lay in between corresponding 
data for p = 0.8 and p = 0.9. Indeed, the most recent estimates read: Pc{p = 0.8) = 
0.2857368(52); p^{p = 0.9) = 0.2492880(35) [i^j, p^{p = 0.8) = 0.2857447(24) [HH]. 

In this study we have calculated the universal ratio F^/F^ for the magnetic 
susceptibility in the critical region 

_ [ F+t"^, T > Tc 

^ " I F-t-T, T < Tc. ^'^^^ 

The singularity of the susceptibility ()24|) is observed for the infinite system only, at 
finite system size L it is rounded-off with the finite maximum value. At each given 
system size L there will be only a finite temperature interval where the susceptibility 
curve overlaps with that for the infinite system. As T approaches the Tf.{oo) the 
finite-size curve deviates from the infinite one. As the result, each estimate for the 
critical amplitudes made at different system sizes will be valid in certain temperature 
range. The other complication for the disordered systems is the presence of some 
distribution of x curves obtained for different disorder realisations (the example for 
the system size L = 48 is presented in fig. IH)). We employed the following method 
for the evaluation of critical amplitudes F+ and F^. The critical temperature for the 
infinite system Tc(oo) was taken as the central one and then the grid of temperatures 
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Figure 6. The example of the distribution for the susceptibihty vs temperature 
curves for different disorder reahsations for the system size L = 48 

around it was considered. The gridness was chosen to be even in scaled temperature 
units V-/''\t\, where t = Kc{oo) — K = ^^^^ff^^^- At each grid temperature the 
separate simulation of 10^ disorder realisations was performed of the length lO^r^; 
and the value for the susceptibility was evaluated. One could, in principle, use the 
histogram reweighting technique to evaluate the intermediate temperatures as well, 
but this will require the study of histograms validity. In this calculation we opted 
for more straightforward approach using the separate simulation runs and haven't 
employed histogram reweighting. 

The values for the critical amplitudes F"*" and T~ can be obtained by plotting the 
scaled susceptibility xl^T- The data are presented in fig. [7|and for the evaluation of 
the ratio /V~ the data points with L^/^\t\ < 1 have been ignored, since they are 
in the FSS regime . This gives the following value for the amplitudes ratio 



0.020 
0.015 
0.010 
0.005 
0.000 ^ 
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r^fi . ! i : ' : 
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Figure 7. The values for the critical amplitudes for the susceptibility and F 
at T > Tc(oo) and T < Tc(oo), respectively 



r+/F- = 1.67 ±0.15 



(25) 
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It compares very well with the recent result F^/F" = 1.62 ± 0.1 [ obtained 
by MC simulations of bond-diluted Ising model at bond dilution ptond = 0.7. The 
theoretical prediction obtained by the field-theoretical renormalization group calcu- 
lations in three-loop approximation reads: F+/F~ = 3.05(32) [EI- As a reason of 
the discrepancy between estimates (j2Sl), [US] and [[T7] one can mention the shortness 
and bad convergence of the series for F"'"/F~, as recognized already in However 
further analysis is required to get a definite answer. In any case, both theoretical 
and MC estimates of the critical amplitudes ratio certainly differ from those of the 
pure 3D Ising model: F+/F^ = 4.70 ^ 4.95 (see [|^ and references therein). 

3. Conclusions and outlook 

In this study we performed computer simulations of the 3D Ising model with 
random site dilution in critical region. Despite the number of simulational studies 
of this model being already performed, some aspects of its critical behaviour and 
simulational details are still awaiting for more detailed analysis. One of these to 
mention is the applicability and the effectiveness of different Monte Carlo algorithms, 
the topic being studied previously for the pure models mainly. However in the case of 
disordered models the impurities may have a profound effect on system clusterisation 
which, in turn, is exploited in cluster methods. 

We examined the most commonly used Metropolis, Swendsen-Wang and Wolff 
algorithms and evaluated both the energy autocorrelation times and the respective 
speed of algorithms in generating statistically independent configurations. We found 
that for the linear lattice sizes up to L = 100 the Wolff algorithm is more preferable 
to use and hence it was applied in our study. 

We used the workstations cluster of the ICMP currently equipped by the Athlon 
MP-2200+ processors. The typical simulation times range from 0.1 seconds per 
1000 MC sweeps for the smallest lattice size L = 10 (Wolff algorithm) to about 953 
seconds for the largest lattice size L = 96 (Metropolis algorithm) per one disorder 
realisation. 

Following the standard ideas of the finite size scaling, the critical indices have 
been calculated. In this way we complemented the existing results. Based on the 
fact that the correction-to-scaling terms appear to be particularly small in the 
region of spin concentration p ~ 0.8, we have performed the simulation for the con- 
centration p = 0.85 only. Obtained values for the critical exponents (fT^. ((201), (EH) 
slightly differ from their asymptotic values as obtained by the field-theoretical RG 
approach [120] and by recent MC simulations [Ej, [HE]- This is the consequence of 
the fact that for the spin concentrations and lattice sizes chosen here the system 
still crosses over to the asymptotic regime. 

Contrary to the critical indices, the values of critical amplitude ratios received 
so far less attention. In our study, we get an estimate for the magnetic susceptibility 
critical amplitude ratio T'^ /r~ = 1.67±0.15. This value is in a good agreement with 
the recent MC analysis of the random-bond Ising model [E] giving further support 
that both random-site and random-bond dilutions lead to the same universality 
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class. 

The comparative analysis of different MC algorithms (including Metropolis, 
Swendsen-Wang and Wolff ones) which is a subject of this study will be used in 
our future work on influence of different types of structural disorder on criticality. 
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